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Abstract 

Node localization plays an important role in many practical applications of wireless underground sen- 
sor networks (WUSNs), such as finding the locations of earthquake epicenters, underground explosions, 
and microseismic events in mines. It is more difficult to obtain the time-difference-of-arrival (TDOA) 
measurements in WUSNs than in terrestrial wireless sensor networks because of the unfavorable channel 
characteristics in the underground environment. The robust Chinese remainder theorem (RCRT) has 
been shown to be an effective tool for solving the phase ambiguity problem and frequency estimation 
problem in wireless sensor networks. In this paper, the RCRT is used to robustly estimate TDOA or range 
difference in WUSNs and therefore improves the ranging accuracy in such networks. After obtaining the 
range difference, distributed source localization algorithms based on a diffusion strategy are proposed 
to decrease the communication cost while satisfying the localization accuracy requirement. Simulation 
results confirm the validity and efficiency of the proposed methods. 

Index Terms 

Wireless underground sensor networks, source localization, Chinese remainder theorem, time-difference- 
of-arrival (TDOA). 

I. Introduction 

Wireless underground sensor networks (WUSNs) are an important extension of terrestrial wireless 
sensor networks, as they can be used to estimate the location of earthquake epicenters, underground 
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explosions, microseismic activities in mines, etc. Normally, the sensor nodes in WUSNs are buried 
underground and they exchange information wirelessly via the dispersive underground channel. Some 
experimental results with WUSNs are reported in Q. 

In terrestrial wireless sensor networks, time difference of arrival (TDOA) measurements are widely 
used for node localization |1I]. Because of the physical characteristics of the dispersive underground 
channel and the heterogeneous network architecture of WUSNs, source localization for WUSNs based 
on TDOA is more challenging 121. In a dispersive medium, we cannot directly obtain the range differences 
between the source and sensors from TDOA measurements since the propagation velocity is a function 
of frequency; different frequency components will have different propagation delays 131. 

Determining the location of sensor nodes is important in many practical applications of wireless 
underground sensor networks. The objective of a positioning system is to determine accurate node 
locations with low complexity and communication cost. Localization algorithms for traditional terrestrial 
wireless sensor networks can be classified into two types: range-based methods and range-free methods. 
Range-based methods usually have higher location accuracy than range-free ones while demanding 
additional hardware cost H. 

In 161, a distributed TDOA estimation method that relies only on radio transceivers without other 
auxiliary measurement equipment was presented. Ultra-wideband (UWB) signaling can be used to accu- 
rately achieve time of arrival (TOA) or TDOA measurements, which has the advantages of low-cost 
and penetrating ability, but also has the weakness of short-range. TDOA based algorithms provide 
high localization accuracy, and represent a practical method for estimating range differences and source 
positions in WUSNs. However, this method also faces many challenges. In particular, limited range 
and directionality constraints decrease the accuracy of range difference estimation. We notice that the 
TDOA can usually be obtained from the measurement of a signal's phase which is susceptible to phase 
ambiguity problems. The Chinese remainder theorem (CRT) offers a closed-form analytical algorithm 
to calculate a dividend from several of its corresponding divisors and remainders, and can be applied 
to solve the ambiguity problem here. In our ranging application, the remainders in the CRT are the 
measured "remainder" wavelengths, the divisors are the measuring wavelengths, and the dividend is the 
range difference to be estimated. However, directly using the CRT is not feasible due to its over-sensitivity 
to noise, i.e., a small error in a remainder can lead to very large error in the estimated dividend. To avoid 
this weakness of the CRT, we propose to use a robust Chinese remainder theorem (RCRT) algorithm to 
estimate the range difference, in which the dividend can still be reconstructed with only a small error 
if the errors on remainders are bounded within a certain level 121. As a result, the range differences or 
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TDOAs can be robustly estimated from noisy measurements in WUSNs by using the RCRT. 

After obtaining the range differences using the RCRT, we can estimate the source position based on 
statistical signal processing methods. For traditional terrestrial wireless sensor networks, TDOA based 
localization algorithms are normally implemented in a centralized way. In the centralized solution, all 
nodes relay their TDOA measurements to a fusion center, which uses a conventional localization algorithm 
to obtain the source position, and then sends the global estimate back to every node. This strategy requires 
a large amount of energy for communications IH and has a potential failure point (the central node). 
Distributed strategies are an attractive alternative, since they are in general more robust, require less 
communications, and allow for parallel processing. To address this limitation of centralized processing, 
we propose distributed source localization algorithms using a diffusion strategy in this paper. Diffusion 
algorithms were proposed in iQl- lfTTl and are applicable to distributed implementations since nodes 
communicate in an isotropic manner with their one-hop neighbor nodes, and no restrictive topology 
constraints are imposed. Thus, the algorithms are easier to implement and are also more robust to node and 
link failures. This approach allows nodes to obtain better estimates than they would without cooperation. 

The rest of the paper is organized as follows. Section II establishes the mathematical model of the 
problem, and derives the proposed ranging method based on the RCRT. Section III gives the distributed 
source localization algorithm based on the diffusion strategy. Simulation results are given in Section IV, 
and conclusions are drawn in Section V. 

II. System model and TDOA estimation via the robust CRT 

Due to the large attenuations in underground environments, experimental results show that underground 
to underground (UG2UG) communication is not feasible at the 2.4GHz frequency band |[T2l . Underground 
communication becomes practical only at lower frequencies. As reported in [|13il and lfT4l . a WUSN system 
operating at 433MHz with a maximum transmit power of +10dBm usually achieves a communication 
range of around one meter for UG2UG communication and more than 30m for underground to above- 
ground (UG2AG) communication. These communication ranges have already exceeded the wavelength of 
the transmitted signal, which results in phase ambiguity when the distance difference is directly calculated 
from the phase. In this section, we first propose a method based on the robust CRT to resolve this phase 
ambiguity when computing the range distance. 

Consider a WUSN with L sensor nodes at known positions (jc/,^,), / - 1, 2, . . . , L, receiving a signal 
from a source at an unknown position {x'',y°) through a dispersive medium, as shown in Fig. [T] The 
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distance between the source and the i-th sensor is 

Ri = ylixi-x-f + iy^-yf. (1) 

The range difference (RD) between the i-th and j'-th receivers, denoted by r,,-, is r,j = 7?, - Rj. 

If the medium is non-dispersive, the propagation delay for any frequency is constant. However, in 
a dispersive medium, the signal propagation velocity is a function of the signal's frequency, denoted 
by Vk for frequency ojk; that is, the propagation delay is frequency dependent. We denote the delay of 
propagation at frequency oj^ for sensor / as 

Tk,i = —■ (2) 

We assume that the source transmits a sinusoid signal Sk(t) = e''^*' at frequency cok- The i-th sensor 
receives the signal as 

skjit) = + mit), (3) 

where nj{t) is the noise at sensor /, and = 1,2, ...,L are independent white Gaussian noise 

processes. Similarly, the received signal at sensor j is 

Sk,j(t) = e''^*''-^*-^' + nj(t). (4) 

By taking the cross-correlation of Skj{t) and Skjit), we get 

C[sk4{t)Tsk,iit)dt 
hu = . (5) 

It is easy to show that 4 ,, is an asymptotically unbiased estimator of JYtJ^ foj- which it follows 

that 

where the range difference r,^ is contained in the phase of 4,,. We denote the phase of Ikjj by 4>k,ij, i-C, 

(f>k,u - mod 27T. (7) 

One can determine r,, from (pkjj- However, there are two issues of concern with this approach: 1) the 
value of (pkji is folded by 2n and 2) the measurements are noisy, i.e., 

(^krij 

= (pkji + 2T:bk Vk, (8) 

Vk 

where bj^ is the quotient (folding integer) and Vk is the noise at frequency ajj^. 
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For issue 1, since 4>k,ij £ [0, 2;r), no matter how large the actual r,y is, the RD r/^,,y = 4'k,ij ■ Vk/cok 
converted directly from (ptjj is always within one wavelength of the signal, which is Ak = Invt/cok- 
Therefore, r/y cannot be determined uniquely from a single (ptjj- We add more constraints to confine the 
solution space by measuring the phase at different frequencies, oj^, k = l,2,...,K. The Chinese remainder 
theorem (CRT) provides a solution to evaluate a dividend from its remainders. We can regard r,; as the 
dividend, Ak as the divisor, and rkjj as the corresponding remainder, i.e., 

rij = bkAk + rkjj. (9) 

CRT tells that a positive integer r,j can be uniquely reconstructed from its remainders rkjj modulo K 
positive integers Ak, if r,; < lcm{Ak},k = l,...,K, where lcm{-) denotes least common multiple. It seems 
that the CRT gives a perfect solution to this problem. However, when we consider measurement noise, 
the traditional CRT |[T6l is not suitable because it is sensitive to noise. A small error in the remainder can 
result in a large error in the estimated dividend. Therefore, we adopt a robust CRT method ||7l. Theorem 
1 in Q proves that the robust CRT can tolerate an error in rkjj that is bounded hy v < B/4, where B is 
the greatest common divisor (GCD) among the divisors. To be more specific, if all the remainder errors 
are not greater than the error bound v, the estimation error of the unknown dividend is upper bounded 
by V, i.e. 

h - hil < V. (10) 

Based on this robust CRT, we provide the following solution to solve the problem: 
Step 1. Estimate (pkjj by calculating the phase of 4 ,,. Denote the estimated phase by 4>k,ii £ [0,2n), and 

the corresponding distance converted from the phase by 

(pk ij 

h,i} = -^h, where r^-,,, e [0, Ak), k = I, K. (1 1) 

ZTT 

Note that the CRT is commonly expressed over the ring of integers while the distance is a real number. 
Therefore, we extend the algorithm from the integers to the reals by introducing a real common factor 
among the divisors as in ifTSl . We choose those communication frequencies so that the Ak,k = 1, ...,K, 
have a real common factor of B which satisfies 

Ak = BFk, (12) 

where Tk are co-prime integers, i.e., (rm,r„) = 1, for 1 < m,n < K,m t n, and (•,•) denotes GCD. 
According to the prior discussion, we have the following equation: 

rij = bkAk + Vk^j ,k= l,...,K, 

(13) 

= bkAk + h,ij + ^rk^j 



December 20, 2011 



DRAFT 



6 



where Ar/; ,, denotes the measurement errors of r^jj with 



< V, and V is the eiTor bound. 



Step 2. For notational convenience, define the following symbols: 



r = ]~|r,, (14) 

k=l 

7k - r/r,. (15) 



Calculate and find the sets, 



Sk = \{bi,bk) = SiTg min \bkAk + hjj - biAi - fiji\\ ,k = 2,3, ...,K, (16) 

where 

n, - ((^1,^^)10 < ^1 < ri - 1,0 < h < 7k - 1} (17) 

is the solution space of the quotients. 

The set can be regarded as the optimal combination of the quotients bi and b^ with which the 
difference of the estimated dividends achieves its minimum. 
Step 3. Let Sk,i denote the first element bi of the 2-tuples, 

Sk.i - {bi\ibi,bk) e Sk for some bk). (18) 

Calculate the intersection set of Sk,i- 



S=f]Sk,i. (19) 



k=2 



In Q, it was proved that if the error bound is less than B/4, the set S contains only the true value of 
bi, i.e., S = {bi}. In addition, bk can also be determined from Sk, that is, if {bi,bk) e Sk, then bk = bk 
for 2 < k < K. Therefore, with all the quotients being determined correctly, the error of the estimated r,y 
is therefore bounded by (fTOl ). where r,j is obtained by averaging the estimates corresponding to different 
wavelengths, i.e., 

1 ^ 

ra = y (bk-ik + hjj). (20) 

Remark: One may notice that to find the set Sk in (fT6l) we need to search among the possible values 
of £lk, which is a 2-D search problem with the order of {r2Tj,...TK)^ and requires a high computational 
complexity. However, the complexity can be decreased by using a fast algorithm proposed in 171 and HI 
to a 1-D search problem with the order of only 2(L - l)r/. One can easily apply the fast algorithm to 
the application in this paper. We do not describe this fast algorithm herein since this is not the focus of 
this paper. 
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III. Diffusion Algorithms for Source Localization 

After obtaining the range differences using the RCRT, a localization algorithm can be used to estimate 
the source position. Herein, we propose an algorithm for source localization that employs diffusion 
strategies proposed in |[9l- lfTll . First, the nodes in a cluster send their measurements to the cluster head. 
Then, the cluster head determines a local estimate using these measurements locally. After obtaining the 
local estimates, cluster heads exchange their local estimates to achieve diffusion. Before describing the 
distributed diffusion process, we first discuss the global WLS estimate when all measurements are sent 
to a fusion center for estimating the source location. 

A. Global Weighted Least Squares Problem 

Consider a set of cluster heads and K = MN sensor nodes (each cluster head is associated with 
M sensor nodes) spatially distributed over some region with known locations x,'s (we consider the 2- 
D Cartesian coordinate system). The objective of the network is to collectively estimate an unknown 
deterministic column vector - the location ;c of a source. All the nodes (cluster heads and sensor 
nodes) in the network can measure the signal transmitted from the source. The sensor nodes transmit 
their measurements to the corresponding cluster head and each cluster head forms a set of M TDOA 
measurements with itself being the reference. Then, cluster heads can send their TDOA measurements to 
the fusion center. The TDOA measurements are formed one at a time by comparing the signal from the 
cluster head and the signal from a sensor node, thus leading to uncorrected estimates if the estimation 
period is longer than the typical coherence time of the mobile radio channel. Finally, the fusion center 
obtains altogether K TDOA based range difference measurements r, /s in the network. 

The scalar model for TDOA based range difference is given by 

r,j = \\x - XiW - \\x - XjW + riij, (21) 
stack all the r,,,'s into a vector and write 

r - s(x) + n 

where t = col{r,_,}, s{x) = col{||.x - Xi\\ - \\x - Xj\\} and n = col{«, Each rij can be alternatively denoted 
as [r]/ to reflect its location in vector r. The corresponding covariance matrix for n is denoted as W and 
lV-diag{(r2,...,(r|}. 
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The global weighted least squares estimator for estimating x given t is thus 

xg = argmin(r - s(x))^W~^(r - s(x)) 



= argmin ^ cr^{[r]i - [s{x)i])\ (22) 



! = 1 



Assuming the «, /s are Gaussian, then the covariance of xq attains the corresponding CRLB, which is 
given by ifTSl 

cov(^g) = (P'^W-^P)-^ (23) 



where 

d[s{x)], 



d[x]j 

and Xq denotes the true position of the source. 

For the TDOA measurements, P is a ^ x 2 matrix and the elements of P are given by 



(24) 



[P]l,l = 

[P]l2 



\\Xi - X\\ \\Xj - x\\ 

[xih - [xh [xjh - [xh 



X=Xo 



\\Xi - X\\ \\Xj - X\\ 

where we assume [^(^c)]/ involves nodes (a sensor node and a cluster head) / and j. 

B. Local Weighted Least Squares Problem 

For each cluster head k, it has access to limited data from its neighbors. It can then solve the WLS 
problem locally as 

K 

Xk = arg min ^ c^crf ([r],- - {s{x)i\f (25) 
1=1 

where c,,i's are the associated weights for node k and c,^, = if [f], is not accessible by node k. Let C 
denote the K x N matrix with elements Cj^. We require l^C = 1^, where 1 denotes the x 1 column 
vector with unit entries. 

A local estimate can also be written as 

Xk - (P'^W-^CkPT^P'^W-^CkZ = LkZ (26) 

where Ck = diag{Ce<.} (e<. is an x 1 vector with a unity entry in position k and zeros elsewhere) and 
z = r - s{x) + Pjc|^=vo- P can be estimated at Xk. 

The covariance matrix associated with Xk is given by coy(xk) = LkWL^. Here Lk contains only local 
information due to the selection ability of Ck- 
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The estimation fusion method can then be used to fuse these local estimates in a distributed manner 
by utilizing the local covariance matrices as proper weights and the estimation fusion method can be 
shown to achieve the performance of the global estimation. However, it involves covariance estimation 
and matrix inversion. 



C. Diffusion Algorithm 

Besides estimator fusion, we can also use a diffusion algorithm to perform distributed estimation. For 
each cluster head k, at the ith time epoch, it exchanges its local estimate with its neighboring cluster 
heads and updates its local estimate using a diffusion algorithm: 



N 

XL 

!=1 



i = ^ ai^khh (27) 



where the a/,i's are the diffusion coefficients. Eq. (|27] ) can be considered as a weighted average of the 
local estimates in the neighborhood of node k. Assume all the local estimates are unbiased. Then in 
order for jci, to be unbiased, we require l^a^ = 1^ where = [ai k, . . . ,ajM^tV . The diffusion process 
is repeated until all the local estimates have converged, i.e., Wx^j+i - x^jW < e Wk, where e is a (small 
positive) design parameter. 

One possible choice for the weights ai^j^ is to consider the degree of connectivity, which is 

deg// 2 deg„, / e Nk 

nm (28) 

0, otherwise 

where deg/ denotes the cardinality of cluster head /'s neighbors (also cluster heads) and A^^. denotes the 

set of neighboring cluster heads of head /. Such choice has been observed to yield good results for the 

diffusion algorithm in general ||9l- 

However, this method does not consider the reliability of different local estimates. The reliability of 

a local estimate is reflected in its associated covariance matrix. But estimating the covariance matrix is 

not an easy task for a nonlinear weighted least-squares estimator (WLSE). Here we propose a method 

for setting appropriate a/^^'s which reflects the reliability of different local estimates to a certain extent 

without requiring use of the covariance matrices. 

Since each local estimate contains certain errors, when sorting them in respective dimensions (each 

dimension of is treated separately), the middle ones are more reliable. Therefore, for cluster head k, 

at the ith time epoch, it first finds the median x^j along respective dimensions among its received local 

estimates xi/s. Then for each obtained local estimate, head k calculates w'^^ = exp(-||.x/, - x^jW^/y) if 
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/ e Nk- Otherwise, wj^. = 0. y is a parameter which controls how rapidly the weight decays as \\xij - x^jW 
grows. The diffusion coefficient aj^ is then set to 

a[k = ^yY.^'jr (29) 

Here we explicitly indicate the time epoch of and aj^ in the superscript. It can be seen that the larger 
the deviation between a local estimate xij and Xj^j, the smaller the weight assigned to this local estimate 
and vice versa. Obviously, l^aj^ = 1^. 

Another method of interest is to use an optimization technique. That is, to set ait^ such that the trace 
of the covariance matrix of x^j is minimized. We start by examining the first time epoch / = 1. Then, 
we have 

N N 

1=1 1=1 

where we have used (l26l ). 

The covariance matrix of x^j is thus given by 

N N 

covixtj) = (J] a[,L!j)WiY, a[,L,jf (31) 
i=\ 1=1 

and its trace is 

N 

tr(cov(^,,,)) = Yj «L,^<,^.ti-(WL^,/^n,) = {a'fQai (32) 

m,n=l 

where a[ - [a\ . . . ,a'j^jj^ and [Q]m,n = ^(WLl^iLnj). To be unbiased, we also require l^a^ = 1^. 

Therefore, to find the optimal in the sense of minimizing the diffusion covariance matrix, we need 
to solve an optimization problem: 

a[ = argmin(4)^e4' ^-f- ^^4 ^ 4 ^ 0' 4,k = V / ^ A^,,. (33) 

This problem is convex if 2 > and thus can be solved fast and efficiently. 

After determining at time epoch /, Lij+i is updated as Lij+i = J^f^^ a'^^Lij and Q will also be updated 
accordingly. Then the above optimization process can be performed iteratively until estimates converge. 

This optimization method can enhance the distributed fusion performance at the cost of slightly higher 
computational complexity compared with simply setting a/i- according to ( [28] ). 

To decrease the computational cost, the optimization process can be implemented only once for the 
first diffusion at each node, and the weight remains the same for the latter diffusions. Since after one 
diffusion, the updated local estimates becomes much less dissimilar and thus the weights will have much 
less influence on the followed diffusions. 
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IV. Simulation results 

In this section, simulation results are given. We first study the ranging performance using the RCRT 
in Subsection A. Then, the localization accuracy is introduced in Subsection B. 

A. Ranging Performance 

Assume that the signals are transmitted at 3 frequencies, i.e., ^"=3, 6=80, ri:={15, 16, 17}, and the corre- 
sponding three dividers are {1200, 1280, 1360} which represent the wavelengths = {120, 128, 136}mm. 

K 

According to the RCRT, the maximum estimation distance is d^ax = ^11 ri.=32640mm. Each trial of 
simulation generates random integers r,j, which are uniformly distributed over [0, JmaxJ- And there are 
1000 trials for each SNR. The result is shown in Fig. jV] 

Then we consider the effect of B on the performance. We set B - 100, 200, 300, respectively and fix 
Tk-{1 ,9]. The result is shown in Fig. 2 with 10000 trials for each SNR. According to the result, changing 
B does not have significant influence on the relative error of distance estimation. 

Next, we compare the performance under different values of F^. with constant B and K. In the simulation, 
we fix B=50, K=3, and let F^. be {7,11, 15}, {29, 33, 37}, {53, 57, 61}, respectively. Fig.3 demonstrates that 
the estimation error increases with F/.. The results from Figs. 2 and 3 can be explained by equations ([TTI ) 
and (fT2l ): the phase measurement error is amplified by B and F^,, i.e., the error of r/. is Ar^ /^ = ^^BTk, 
where A0<._,y denotes the phase measurement error of ^i^^j. The larger Fi- is, the larger error results. 
However, B does not affect the performance because the robustness of the algorithm also linearly increases 
with B which cancels the performance deterioration of the increased B. (The error tolerance of the 
algorithm is < B/4) 

In addition, we consider the scenario in which different sets of F^- are compared under the constraint 

K 

of constant maximum range dmax = B Y\^k- In Fig- 4, We choose Fi.={7, 11, 15}, {5, 11,21}, {3, 11,35}, 

/t=i 

respectively. The simulation results suggest that the performance is better if the differences between the 
Fi- are smaller. 

Fig. 5 demonstrates that using more wavelengths results in better ranging performance. In this simu- 
lation, we fix the maximum estimation distance <imax> and vary K. We consider the cases when K=2,3,4, 
respectively, with B=50 and t/max=3465mm. Fit are set to {33,35}, {7, 11, 15}, {3,5,7, 11}, respectively. 
Simulation results demonstrate that our RCRT based ranging scheme can estimate the range differences 
with high accuracy. 
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B. Comparison Results for Distributed Source Localization 

We now show simulation results for source localization. The algorithms compared are: 1) Global: the 
global weighted least square estimator (l22l ). 2) Diff (con): the diffusion algorithm with coefficients set 
by considering connectivity ( |28] ). 3) Diff(wei): the diffusion algorithm with coefficients set by weighting 
( [291 ). 4) Diff (opt): the diffusion algorithm with coefficients set by optimization (1331 ) and 5) Local: simple 



average of aU the local estimates, i.e., x/i/N. 

The root mean square error (RMSE) is used as a performance metric, which is defined as yjE(\\x - 
The Cramer-Rao Lower Bound (CRLB) on the RMSE based on the entire data (equals to Vtr(cov(I^ 
for Gaussian measurement noise) is also presented as a benchmark. Each simulated point is averaged 
over 200 runs. 

The simulated network consists of cluster heads that are regularly deployed at grid points. The 
distance between neighboring cluster heads is set to 50 (the units are meters, here and below). Each 
cluster head has M associated sensor nodes which are distributed uniformly around the corresponding 
cluster head. The source location is fixed at [60, 70] in all the simulations. The TDOA measurements are 
generated according to ([2T[ ) with ?i,_,'s being Gaussian noises. W is set to W = cr^I with cr = 1. y is set 
to 1. The initial point for the WLSE is always set as the center of the deployment area. 

Here we consider the scenario in which each cluster head exchanges its local measurements with its 
neighboring cluster heads to perform a local estimate, and then exchanges its local estimate with its 
neighboring cluster heads to perform diffusion until convergence. C is set to 

ci,k, [?]/ involves cluster head 1,1 e Nk and I k 

Ck^k, [?]/ involves cluster head k (34) 

0, otherwise 



where 



Clk 



l/max{deg,,degj.}, / e A^^, I k 
1 - ^ / ^ k 



(35) 



0, otherwise. 
Here, ci^k denotes the weight assigned with respect to cluster heads / and k, and c, represents the weight 
assigned to the /th measurement used by the cluster head k. 

First, we fix = 16 and vary the value of M, which changes from 5 to 50 with a step size of 5. The 
RMSEs of respective algorithms are shown in Fig. [7J It can be observed that Global is the best which 
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can attain the CRLB, local is the worst and Diff (con) is always better than local. In general, Diff (opt) 
is better than Diff (wei), and Diff (wei) is better than Diff (con). The performance improvement of Diff 
(opt) and Diff (wei) compared with Diff (con) comes from the consideration of the reliability of the local 
estimates. Though the diffusion algorithms are always worse than Global, the performance differences 
are not significant. As M grows, all the algorithms perform better. 

Fig. [8] shows the corresponding average CPU times of respective algorithms except Diff (opt) whose 
CPU time is typically 10 times that of Global due to its numerical optimization nature (the same 
below). It can be seen that Diff (con) is very time efficient with a CPU time almost the same as 
local. Diff (wei) consumes a little more CPU time than Diff (con), while the CPU time of Global is 
much larger. This demonstrates the advantage of the diffusion algorithm in terms of CPU time and 
computational complexity. Furthermore, we can say the diffusion algorithm has lower communication 
cost and computation complexity than the centralized solution, while the localization accuracy is close 
to the centralized method. 

Fig. |9] shows the average number of iterations before convergence of the three diffusion algorithms. 
We can observe that Diff (opt) requires many fewer iterations compared with Diff (con) and Diff (wei). 
Diff (wei) needs slightly more iterations than Diff (con), which may explain our observation of a slightly 
longer CPU time consumed by Diff (wei). 

We now fix M = 10 and examine the effects of the number of cluster heads on the performance 
of the algorithms. We vary from 4 to 36. The RMSEs of the algorithms are shown in Fig. [TOl It 
is clear that adding more clusters (and thus sensor nodes) will not necessarily improve the estimation 
performance. This is known as the geometric effect of the localization problem. Since the source location 
is outside the convex hull formed by respective added clusters, the corresponding local estimates are not 
good enough. Diff ( con ) and Diff ( wei) thus also show performance degradation. However, Diff ( wei) is 
much better than Diff ( con ) when more clusters are added while the sensor nodes associated with each 
cluster head is fixed. The performance of Global is almost unchanged as becomes large. Due to the 
consideration of the estimation covariance matrices, Diff (opt) shows very good performance. However, 
it will consume much more CPU time. The corresponding average CPU times and average number of 
iterations are shown in Figs. [TT] and [T2l respectively. Similarly, Diff (con) is the most time efficient and 
Diff (opt) requires the smallest number of iterations. 

Then we examine the effect of a on the performance of the algorithms. We set = 16 and M = 10. The 
RMSEs of respective algorithms are shown in Fig. [13] The relative performance of respective algorithms 
is the same as before. As cr enlarges, all of them show performance degradation. The average CPU times 
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and average number of iterations among different algorithms have the same relationship as before and 
thus the figures are not shown here. 

Finally, we examine the choice of y on the performance of Diff (wei). We setA^=16, M^IO and 
cr = I. We generate 200 realizations of the overall TDOA measurement vector, then store and use them 
for all the corresponding simulations. The corresponding RMSEs are shown in Fig. [141 It can be seen 
that an optimal y exists which results in the minimum RMSE for Diff (wei). However, in a large range 
of the choice of y, Diff (wei) can generate better results than Diff (con). 

V. Conclusion 

We have presented energy efficient localization schemes that can achieve high locahzation accuracy 
in wireless underground sensor networks. These distributed localization algorithms require low compu- 
tational complexity and energy consumption based on a diffusion strategy. An accurate RCRT based 
ranging scheme using TDOA to determine range differences between sensors and source that does not 
require time synchronization is also proposed. It has been shown via simulation results that the proposed 
localization algorithms achieve excellent localization accuracy with lower communication cost. In future 
work, we plan to implement our localization scheme in a testbed and verify its performance with an 
actual WUSN. 
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Fig. 1. System model for WUSN localization. 
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Fig. 2. Relative estimation eiTor decreases with increasing SNR. 
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Fig. 3. Comparison between the relative estimation errors for different values of B. 
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Fig. 4. Comparison of the relative estimation error for different values of Ft, when B and K are fixed. 
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Fig. 6. Comparison of the relative estimation error for different values of K. 




Fig. 7. RMSE versus M when A' = 16. 
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Fig. 8. CPU time versus M wlien A' = 16. 




Fig. 9. Average number of iterations before convergence versus M when N = 16. 




Fig. 10. RMSE versus A' when M = 10. 
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Fig. 11. CPU time versus N when M = 10. 
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Fig. 12. Average number of iterations before convergence versus A' when when M = 10. 




Fig. 13. RMSE versus cr. 
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